##########################################################################################

library(data.table)
library(optparse)
library(bedtoolsr)
library("MutationTimeR")

##########################################################################################

option_list <- list(
    make_option(c("--sample"), type = "character"),
    make_option(c("--ccf_file"), type = "character"),
    make_option(c("--vcf_file"), type = "character"),
    make_option(c("--cnv_file"), type = "character"),
    make_option(c("--purity_file"), type = "character"),
    make_option(c("--chr_len"), type = "character"),
    make_option(c("--out_path"), type = "character"),
    make_option(c("--bedtools_path"), type = "character")
)

if(1!=1){
    
    sample <- "JZ423N"
    vcf_file <- "/public/home/xxf2019/20220915_gastric_multiple/dna_combine/mutationTime/input/JZ423N.vcf"
    cnv_file <- "/public/home/xxf2019/20220915_gastric_multiple/dna_combine/mutationTime/input/JZ423N.seg"
    purity_file <- "/public/home/xxf2019/20220915_gastric_multiple/dna_combine/titan/Purity_titan.final.tsv"
    out_path <- "/public/home/xxf2019/20220915_gastric_multiple/dna_combine/mutationTime/result"
    chr_len <- "/public/home/xxf2019/ref/seq/GRCh37.fa.fai"
    ccf_file <- "/public/home/xxf2019/20220915_gastric_multiple/dna_combine/titan/chat/JZ423N_CHAT.txt"
    bedtools_path <- "/public/home/xxf2019/tools/StandTools"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample <- opt$sample
vcf_file <- opt$vcf_file
cnv_file <- opt$cnv_file
purity_file <- opt$purity_file
ccf_file <- opt$ccf_file
chr_len <- opt$chr_len
out_path <- opt$out_path
bedtools_path <- opt$bedtools_path

options(bedtools.path = bedtools_path)


###########################################################################################
## 构建38版本的基因组长度
dat_chrLen <- data.frame(fread(chr_len))
dat_chrLen <- dat_chrLen[1:23,]

dat_chrLen  <- data.frame( seqnames = dat_chrLen[,"V1"] , start = 1 , end = dat_chrLen[,'V2'] , 
    width = dat_chrLen[,'V2'] - 1 , strand = "*" )

chrOffset <- c(0) 
for(i in 2:dim(dat_chrLen)[1]){
    tmp_len <- dat_chrLen[(i-1),"width"] + 1
    chrOffset <- c(chrOffset , chrOffset[i-1] + tmp_len)
}
names(chrOffset) <- dat_chrLen$seqnames
regions <- GRanges(dat_chrLen)
refLengths <- regions
## MutationTimeR:::refLengths[1:23]

###########################################################################################
###########################################################################################

## VCF
## Point mutations, needs `geno` entries `AD` and `DP` or `info` columns t_alt_count t_ref_count.
vcf <- readVcf(vcf_file)

## Purity
dat_purity <- fread(purity_file)
purity <- subset(dat_purity , Sample == sample)$Purity
###########################################################################################

###########################################################################################
## CNV
## Copy number segments, needs columns  major_cn, minor_cn and clonal_frequency of each segment
## clonal_frequency=purity
## https://github.com/gerstung-lab/MutationTimeR/issues/5
## 对于克隆拷贝数片段，clonal_frequency =纯度。在本例中，所有片段都是克隆的。
## 如果存在亚克隆拷贝数片段(如。一个在GRanges对象中有两个坐标相同的条目，但是major_cn和/或minor_cn以及clonal_frequency值不同。
## 该区域的两个clonal_frequency值的和需要等于纯度。
## 请注意，从先验上看，尚不清楚这两种构型中哪一种是祖先状态，哪一种是派生状态(即最近的状态)。
## 通常我们会很自然地假设二倍体是祖先的状态，但并不一定要像WGD的情况那样。对于点突变，通常不指定祖先状态，因为它被隐式假定为匹配的正常状态。
## https://github.com/gerstung-lab/MutationTimeR/issues/5

seg <- data.frame(fread(cnv_file))
seg  <- data.frame( seqnames = seg[,"Chromosome"] , start = seg[,'Start'] , end = seg[,'End'] , width = (seg[,'End'] - seg[,'Start']) ,
                strand = "*" , 
                major_cn = seg[,'Corrected_MajorCN'] , minor_cn = seg[,'Corrected_MinorCN'] , total_cn = seg[,'Corrected_MajorCN'] + seg[,'Corrected_MinorCN'] ,
                clonal_frequency = purity)
seg <- subset(seg , minor_cn!="NA")

## CNV未覆盖的区域改为二倍体，以保证所有的突变都能纳入（若无CNV覆盖，软件不会判断计算该突变的CCF）
## bedtools取未覆盖的区域，起始位置 + 1 ，终止位置 -1 ，保证不会和原来的Seg有重叠
seg_nocnv <- bt.subtract(a = dat_chrLen , b = seg )
dat_seg_nocnv <- data.frame( seqnames = seg_nocnv[,"V1"] , start = seg_nocnv[,'V2'] + 1 , end = seg_nocnv[,'V3'] -1, 
        width = (seg_nocnv[,'V3']  - seg_nocnv[,'V2'] - 2) ,
        strand = "*" , 
        major_cn = 1 , minor_cn = 1 , total_cn = 2 ,
        clonal_frequency = purity)

## 排除长度为1的区域
dat_seg_nocnv <- subset( dat_seg_nocnv , width > 1)

seg <- rbind(seg , dat_seg_nocnv)
## 去除性染色体
seg <-  subset(seg , seqnames!="X")

seg <- GRanges(seg)

###########################################################################################

###########################################################################################
## clusters
## https://github.com/gerstung-lab/MutationTimeR/issues/6
## In the absence of detailed subclonal cluster information, it may be worth introducing a placeholde for subclonal mutations at 50% purity using
# clone_num <- round(0.9 * dim(vcf)[1])
# subclone_num <- round(0.1 * dim(vcf)[1])
## 可估计纯度样本的设置
## clusters <- data.frame(cluster=1:2, 
##                        proportion=max(seg$clonal_frequency,na.rm=TRUE) * c(1, 0.5), ## enabling subclonal mutations at approx. 50% purity
##                        n_ssms=c( 90 , 10)) ## assuming a prior expectation of 10% subclonal mutationss

## 无CNV改变的样本，
## 早期/晚期无法判断
## 观察vaf，基本都不到0.1，但默认cluster判断90%均为克隆突变，不合理，大部分应为亚克隆突变
## 假设90%的为亚克隆
clusters <- data.frame(cluster=1:2, 
                        proportion=max(seg$clonal_frequency,na.rm=TRUE) * c(1, 0.5), ## enabling subclonal mutations at approx. 50% purity
                        n_ssms=c( 10 , 90)) ## assuming a prior expectation of 90% subclonal mutationss


###########################################################################################

###########################################################################################

mt <- mutationTime(vcf, seg, clusters=clusters)

###########################################################################################
## These probabilities are the basis for the following simple clonal states
# table(mt$V$CLS)
vcf <- addMutTime(vcf, mt$V)

## Timings of copy number gains
## The relevant columns are time with 95% confidence intervals time.lo and time.hi as well as the counterparts time2/time2.lo/time2.hi for the 
## second gain in cases where one allele has more than one gained copy. 
## The field time.star indicates a tiers: *** indicates gains +1. ** gains +2, which are found to be slightly less reliable and need certain assumptions about their temporal sequence. * are subclonal gains which are hit an miss.

## 相关的列是带有95%置信区间的time.lo and time.hi，还有对应的time2/time2.lo/time2.hi表示第二个增益，在一个等位基因有一个以上获得的拷贝的情况下。
## time.star表示等级:***表示增益+1，**增益+2（这被发现是稍微不可靠的，需要对其时间序列的某些假设）。
## *是无法定义的次克隆增益。

## head(mt$T)
mcols(seg) <- cbind(mcols(seg),mt$T)


###########################################################################################
#3 info(header(vcf)) <- rbind(info(header(vcf)),MutationTimeR:::mtHeader())
tmp <- data.frame(info(vcf))

dat_ccf <- data.frame(fread(ccf_file))
rownames(dat_ccf) <- paste0(dat_ccf$Chr , ":" , dat_ccf$Start_Position , "_" , dat_ccf$REF , "/" , dat_ccf$ALT )

dat_ccf <- merge(dat_ccf , tmp , by="row.names")
image_name <- paste0(out_path , "/" , sample , "_CCF_mutTime.tsv")
write.table(dat_ccf[,-1] , image_name , sep = "\t" , row.names = F , quote = F )

image_name <- paste0(out_path , "/" , sample , "_CNV_mutTime.seg")
out_seg <- data.frame(seg)[,c(1:3,6:9,11:17)]
# out_seg <- data.frame(seg)
write.table(out_seg, image_name , sep = "\t" , row.names = F , quote = T )

image_name <- paste0(out_path , "/" , sample , "_CNV_mutTime.rda")
save( out_seg , file = image_name)

###########################################################################################
## 画图
image_name <- paste0(out_path , "/" , sample , "_mutTime.pdf")
pdf(image_name , width = 10)
plotSample(vcf , seg , regions = regions )
dev.off()


## This shows the observed and expected variant allele frequencies of point mutations on the top. 
## This is very useful to spot inconsistencies with purity and copy number configuration. 

## As a rule of thumb the states (horizontal bars) should run right through the middle of the clouds of point mutations. 
## Colours indicate the timing category: Blue = clonal [other], purple = clonal [late], green = clonal [early], red = subclonal.

## The middle plot shows the copy number as stacked barplots. Subclonal CN is indicated by fractional bars. Dark grey is major, light grey minor allele.
## The bottom plot shows the estimated mutation time of primary and secondary gains (shaded). Boxes denote 95% CIs. 
## The histogram at the right shows the distribution of timing events. 
## Blue = mono-allelic gains (N:1), pink = CN-LOH/gain+loss (N:0) and green = bi-allelic gains (N:2).

## 这显示了顶端点突变的观察和预期的变异等位基因频率。这对于发现纯度和拷贝数配置的不一致是非常有用的。
## 根据经验，状态(水平条)应该正好穿过点突变云的中间。颜色表示时间类别:蓝色=克隆[其他]，紫色=克隆[晚期]，绿色=克隆[早期]，红色=亚克隆。
## 中间的图显示复制数为堆叠的barplot。亚克隆CN用分数条表示。深灰色为主要等位基因，浅灰色为次要等位基因。
## 下图显示了估计的一次和二次增益的突变时间(阴影部分)。方框表示95%的CIs。
## 右边的直方图显示了计时事件的分布情况。蓝色=单等位基因增益(N:1)，粉色= CN-LOH/增益+损失(N:0)，绿色=双等位基因增益(N:2)。




## 默认所有拷贝都是克隆
if(1!=1){
    ## 若拷贝数未改变，认为是所有细胞都有的主克隆
    seg[which(is.na(seg$clonal_frequency)),"clonal_frequency"] <- 1

    ## 亚克隆片段拆成两个,一个为突变的CNV一个为正常的CNV,保证两个的clonal_frequency之和为1
    add_seg <- c()
    for(j in 1:dim(seg)[1]){
        if(seg[j,"clonal_frequency"]==1){
            seg[j,"clonal_frequency"] = purity
        }else{
            ## CNV改变
            seg[j,"clonal_frequency"] = round(purity * seg[j,"clonal_frequency"] , 2 )
            ## CNV不变
            tmp <- data.frame( seqnames = seg[j,"seqnames"] , start = seg[j,'start'] , end = seg[j,'end'] , 
                width = (seg[j,'end'] - seg[j,'start']) ,
                strand = "*" , major_cn = 1 , minor_cn = 1 , total_cn = 2 ,
                clonal_frequency = purity - seg[j,"clonal_frequency"]
             )
            add_seg <- rbind(add_seg,tmp)
        }
    }

    seg_out <- rbind(seg , add_seg)
    seg <- seg_out
}
